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The boreal forest accounts for one-third of global forests, but remains largely inaccessible to ground-based 
measurements and monitoring. It contains large quantities of carbon in its vegetation and soils, and research 
suggests that it will be subject to increasingly severe climate-driven disturbance. We employ a suite of 
ground-, airborne- and space-based measurement techniques to derive the first satellite LiDAR-based esti- 
mates of aboveground carbon for the entire circumboreal forest biome. Incorporating these inventory tech- 
niques with uncertainty analysis, we estimate total aboveground carbon of 38 ± 3.1 Pg. This boreal forest 
carbon is mostly concentrated from 50 to 55° N in eastern Canada and from 55 to 60° N in eastern Eurasia. 
Both of these regions are expected to warm >3 °C by 2100, and monitoring the effects of warming on 
these stocks is important to understanding its future carbon balance. Our maps establish a baseline for future 
quantification of circumboreal carbon and the described technique should provide a robust method for future 
monitoring of the spatial and temporal changes of the aboveground carbon content. 

© 2013 The Authors. Published by Elsevier Inc. All rights reserved. 


1. Introduction 

The boreal forest biome covers 12.5 ±1.5 million square kilometers 
(km 2 ) (Dixon et al„ 1994; Wulder, Campbell, White, Flannigan, & 
Campbell, 2007) accounting for -27% of the global forested area, and 
stores a vast quantity of carbon in its vegetation and soils (Pan et al„ 
201 1 ). This biome has experienced the greatest magnitude of warming 
of all biomes over recent decades with a 1-3 °C, or more, increase in 
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seasonal surface air temperature (Hansen et al., 2006), and has experi- 
enced the greatest proportion of forest cover lost of all forested biomes 
from 2000 to 2005 (Hansen, Stehman, & Potapov, 2010). There are large 
uncertainties in estimates of the size and spatial distribution of boreal 
forest carbon (C) stocks, and we do not know how it will respond to fu- 
ture warming, forcings and feedbacks (Grosse et al., 2011; Tchebakova, 
Parfenova, & Soja, 2011). Underlying this biome is rich organic soil C 
containing 50% of the total global soil pool (1642 Pg C) that has accumu- 
lated over the past 10,000 + years following glacial retreat, with 88% of it 
locked in a permafrost C pool (Tarnocai et al„ 2009). Warming-induced 
thawing of permafrost could enhance decomposition and increase C 
flux to the atmosphere, negating the long-term sink this biome has sup- 
ported for centuries (Canadell et al., 2007). To date, the spatial distribu- 
tion of aboveground C stock throughout this biome has not been 
characterized with satellite LiDAR (Light Detection and Ranging) mea- 
surements. Furthermore, there is no universally recognized definition of 
the boreal forest, and this can lead to errors in C stock estimates 
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(Houghton et al., 2007). To reduce uncertainties in the global C balance, 
we must improve our characterization of C stocks in the boreal forest. 
This is key to improving our understanding of how this biome may 
change with warming climate. 

The boreal forest biome extends in a circumpolar band 13,400 km 
in length around the Northern Hemisphere roughly between 45° N 
and 70° N. It is bound by tundra to the north and by temperate decid- 
uous forests or savanna/prairie/steppes to the south. Coniferous 
species of spruce, pine, and fir as well as deciduous larches, birches, 
alders, and aspens dominate vegetation cover. Forest fires and insect 
outbreaks are the dominant agents of disturbance and the size of 
these events dwarfs those in other biomes because they are typically 
not suppressed as most do not impact population centers (Stocks et 
al., 2001 ). Harvest occurs primarily in southern portions of the boreal 
forest due to higher productivity and limited access in the north. All of 
these disturbances have had a major impact on C throughout Canada 
(Kurz & Apps, 1999) with similar impacts during the 1990s in Siberia 
(Kharuk, Ranson, & Dvinskaya, 2008). Recently in western North 
America, mountain pine beetle outbreaks ( Dendroctonus ponderosae ), 
facilitated by fire suppression and changes in climate, have resulted in 
widespread tree mortality that has markedly reduced C stocks pres- 
ent in the impacted forests (Kurz et al., 2008). In eastern Canada, out- 
breaks of spruce budworm (Choristonneura fumiferana Clem.) and 
hemlock looper (Lcimbdina fiscellarici Guen.) could potentially negate 
the forest C sink (Dymond et al., 2010). The most recent catastrophic 
insect outbreak in northern Eurasia occurred in the mid 1990s 
(Kharuk, Ranson, Kozuhovskaya, Kondakov, & Pestunov, 2004). Cli- 
mate warming may have already induced more fire disturbance and 
frequent insect outbreaks extending the zone northward, and may 
have even begun a biome shift into the taiga-tundra transition zone 
(Ropars & Boudreau, 2012). Increased productivity (Goetz, Mack, 
Gurney, Randerson, & Houghton, 2007) and increased woody vegeta- 
tion cover containing more C have already been observed in the high 
arctic tundra (Ropars & Boudreau, 2012). The vast area of the boreal 
forest, its large C stocks, and the rapid change that is occurring, all 
provide a rationale for production of a baseline map of C stocks 
throughout the biome. Here we report aboveground C total mass 
and density for 2005, with estimates of uncertainty, developed 
using model-based estimators (Stahl et al., 201 1 ), by leveraging forest 
structure data obtained by the Geoscience Laser Altimeter System 
(GLAS) onboard the Ice, Cloud, and land Elevation Satellite (ICESat). 
An underlying assumption of these estimators is that the regression 
models that relate ground-based estimates of biomass to airborne 
LiDAR or satellite LiDAR measurements must be correctly specified 
for the area of interest, otherwise estimates may be biased. In addi- 
tion, descriptions of procedures which may be used to measure re- 
gional C stocks in remote and inaccessible regions of the boreal 
forest where little information exists are provided. 

2. Data and methods 

Spaceborne LiDAR remote sensing provides an opportunity to con- 
duct forest C assessments because of the capability for repeated mea- 
surements over large areas of interest. The launch of the GLAS sensor 
onboard ICESat in January 2003 offered the first opportunity to eval- 
uate the potential of a spaceborne LiDAR for measuring structural 
properties of the boreal forest. ICESat stopped collecting data in 
2009, although it has provided a wealth of data that was used to pro- 
duce global forest height maps (Lefsky, 2010; Los et al„ 2012; Simard, 
Pinto, Fisher, & Baccini, 201 1 ). GLAS estimates have also been used for 
large area stand volume and/or C in the tropics (Baccini et al., 2012; 
Saatchi et al., 2011), North America (Nelson, Boudreau, et al., 2009), 
and Siberia (Nelson, Ranson, et al., 2009). Other studies have related 
satellite-derived normalized difference vegetation index (NDVI) 
data (Dong et al„ 2003) and spaceborne synthetic aperture radar 
(SAR) (Santoro et al., 2011) to ground-based estimates of biomass. 


Our study used a hierarchical sampling approach to estimate treed 
circumboreal forest carbon (C); specifically from total aboveground 
live dry biomass density (AGB); we approximate C to be 50% of AGB 
(Houghton et al., 2000). This hierarchical scheme related AGB data 
collected on ground plots to airborne and spaceborne LiDAR measure- 
ments. The goal was to develop model-based estimates of boreal for- 
est AGB for land cover strata within ecoregions for five circumpolar 
regions — Alaska, western Canada, eastern Canada, western Eurasia, 
and eastern Eurasia. We utilized two similar approaches in North 
America and Eurasia. We based these methods on available data and 
lessons learned from prior investigations using GLAS-based models 
to predict AGB and C stock. These data have been used in two smaller 
regional studies in 2004 for a portion of Siberia (Ranson, Sun, Kovacs, 
& Kharuk, 2004a, b), in 2005-2006 for southern Norway (Nassset et 
al., 2011; Stahl et al., 2011), and in 2005 for the province of Quebec 
in Canada (Boudreau et al., 2008; Nelson, Boudreau, et al., 2009). 
We supplemented existing ground and airborne data collected in 
prior studies, with new field measurements and additional airborne 
LiDAR data. 


2.1. Ground plot data 

The Canadian Forest Service (CFS) provided access to relevant 
plot data holdings and liaised with provincial and territorial forest 
resource management agencies to obtain access to geo-located 
ground plots in Canadian boreal ecoregions. Based on the large area 
and number of jurisdictions involved, numerous individual re- 
searchers, provincial foresters, and industrial foresters measured 
these plots. Field measurements used in this study were collected 
in Northwest Territories (2006-2008), Saskatchewan (2004-2006), 
Ontario (2006-2007), and Quebec (2001-2004). The CFS has 
established species-specific, national-level equations (Lambert, 
Ling, & Raulier, 2005) that we used in this study to convert ground 
plot measurements to AGB. 

In Alaska, we accessed 361 geo-located ground plots measured by 
the following U.S. organizations: Forest Service — Kenai Peninsula; De- 
partment of Defense — military installations near Fairbanks; and the 
National Park Service — near Denali N.P., Wrangell — St. Elias N.P., and 
Yukon — Charley Rivers National Preserve from 2005 to 2007. United 
States Forest Sen/ice (USFS) in conjunction with a related study com- 
piled the plot locations and associated ground measurements. 

Field data were collected with a variety of techniques across Eurasia 
in 2007, 2008 and 2010. In western Eurasia 201 ground plots were col- 
lected in a 960 km 2 area in south-eastern Norway. More information 
about how these data were collected can be found in Ntesset et al. 
(2011). In southern Siberia and northeast China, 322 fixed area ground 
plots were established on GLAS pulse centroids using a Trimble GeoXT 
differential global position system (dGPS). Of these 322 plots, we ulti- 
mately used 55 to establish the eastern Eurasian model relating ground 
estimates of biomass to GLAS measurements. We discarded a majority 
of the unused ground plots either because the associated GLAS mea- 
surements were outside of our acceptable temporal window or because 
of insufficient signal. The 55 ground plots allowed us to develop predic- 
tive models that related ground-measured biomass to GLAS measure- 
ments. Before field data collection, GLAS waveforms were visually 
examined to select for a strong vegetation signature for potential field 
measurement, not confounded by clouds or slope. We selected GLAS 
shots for which data suggested a strong vertical signature of vegetation. 
In northern Siberia, sample plots were collected along the Kochechum 
River in 2007 and the Kotuykan River in 2008. Both expeditions includ- 
ed transport into the field by helicopter, then boat transport down-river 
accessing GLAS shots within 2-3 km of the river for -2-week long expe- 
ditions. During the summer of 201 0, GLAS shot locations were accessed 
by four-wheel drive vehicle near Chylum-Ket River region in the west- 
ern plains of Siberia. 
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We adopted a two-step procedure to calculate biomass at GLAS 
shot centroids. Field measurements and observations, made in single 
radius (circular) plots of 10 or 15 m consisted of: 1) tree species and 
diameter at breast height (DBH), ±0.1 cm, for all trees > 3.0 cm in 
the entire plot; and 2) a sampling of tree height measurements 
from small, medium, and tall trees for each plot to characterize the 
range of heights. In general, a 10 m radius plot was employed in 
dense southern stands (between 50° and 60° N) and a 15 m radius 
plot was employed in sparser stands above 60° N latitude. We were 
able to collect field data across a broad area of Siberia; but we ac- 
knowledge a notable geographic gap in our Eastern Hemisphere 
ground measurements. Access is limited to these forests due to cost 
and available field personnel to collect measurements over vast re- 
mote areas, particularly in the Russian Federation Far East. 

2.2. Airborne LiDAR data 

2.2.1. PALS 

The Portable Airborne Laser System (PALS; Nelson, Parker, & Horn, 
2003) instrument serves as an intermediate sampling tool in North 
America to extend the spatially limited Canadian National Forest In- 
ventory (NFI) and the U.S. Forest Inventory Analysis (FIA) ground 
measurements to the continental-scale GLAS observations. The sys- 
tem records the LiDAR range/amplitude data stream interleaved 
with dGPS location information, and also records a dGPS-annotated 
video of targets overflown (Nelson et al., 2003). The airborne system 
has been used to measure forest height and canopy closure along 
GLAS orbital transects 1000s of km in length. The nominal flight pro- 
file for the Alaska and Canada airborne laser data collection required 
the pilots to fly -150 m-200 m above-ground-level (AGL) at a 
speed of approximately 97 kn (-180 km h 1 or 50 m s~'). Eastern 
Canada data (Quebec) were collected in the summer 2005, central 
(Ontario) and western Canada (Northwest Territories and Saskatche- 
wan) were collected in the summer of 2009, and Alaska data were 
collected in the summer of 2008. Data were typically acquired at 
333 Hz, yielding an effective along-track post spacing of 15 cm. The 
size of the laser spot at target was -45 cm, so there was (intentional- 
ly) significant oversampling along track. Post-flight processing was 
done to identify ground returns. A spline was fit to the ground points 
to define a ground line. Once the ground line was identified, canopy 
height could be calculated for every pulse. The PALS data stream con- 
tains ranging information that can be used to calculate canopy height, 
top-of-canopy height variation, and canopy closure (see Nelson et al., 
2003, Fig. 1). 

2.2.2. ALTM 

In Southeast Norway airborne data were collected with an Optech 
Airborne Laser Terrain Mapper (ALTM) 3100 system flying at an alti- 
tude of -1850 m at a speed of approximately 145 kn (-270 km h" 1 
or 75 m s~’). Data were acquired in the summer of 2005 with a 
pulse repetition frequency of 50 kHz, with a scan frequency of 
71 Hz, resulting in a point density on the ground of approximately 
0.7 m 2 . The maximum scan angle was 15° but pulses emitted at an 
angle >13° were discarded during subsequent data processing. 
These airborne data served as an intermediate sampling tool to ex- 
tend spatially limited survey plots. More information about how the 
data were processed can be found in Naesset et al. (2011). 

2.3. ICESat-GLAS LiDAR data 

The Ice, Cloud, and land Elevation Satellite (ICESaT) was launched 
in January 2003 with the GLAS instrument on-board. We accessed 
and processed GLAS GLA-01 (waveform) and GLA-14 (land) data dis- 
tributed by the National Aeronautical and Space Administration 
(NASA) Goddard Space Flight Center (GSFC) (http://reverb.echo. 
nasa.gov). We used the satellite data as a sampling tool to understand 


Table 1 

ICESAT-GLAS periods of data used in this study listed by laser acquisition. 


Laser 

identifier 

Start date 
(mm-dd-yr) 

End date 
(mm-dd-yr) 

Pulse 

footprint size 

L2a 

09-24-03 

10-15-03 

95 x 52 m 

L3a 

10-03-04 

10-15-04 

61 x 47 m 

L3c 

06-08-05 

06-13-05 

61 x 47 m 

L3f 

06-08-06 

06-26-06 

61 x 47 m 


forest vertical structure (Harding & Carabajal, 2005; Sun, Ranson, 
Kimes, Blair, & Kovacs, 2008). GLAS carried three lasers that were 
used sequentially to collect measurements on 33-55 day sub-cycles 
from January 2003 through October 2009. In our study, we used 
GLAS pulses from L2a, L3a, L3c, and L3f from the early fall of 2003 
through late spring of 2006 (Table 1). 

Cycles of data were acquired in June and September-October. Our 
study focused on within-growing-season data, and when available 
leaf-on conditions. The lasers of GLAS had a slightly elliptical ground 
pulse spots occurring every 172 m along the ground track of the 
spacecraft (Abshire et al., 2005). GLAS is a waveform instrument 
(Harding & Carabajal, 2005), unlike PALS, which records a first-last 
return range measurement from aircraft to target for each pulse. 
GLAS records the brightness of the 1.064 pm, near-infrared return in 
one-nanosecond increments as the pulse traverses from the top of 
the target to the ground. Over trees, the sequential returns recorded 
for a single pulse provide an initial return from the top of the canopy. 
Through sequential secondary returns in 15 cm vertical bins, they 
also provide ranging measurements to sub-canopy layers and the 
ground as the pulse traverses vertically from top to bottom (Ranson 
et al., 2004a). 

Each individual waveform can be analyzed to extract a number of 
measurements related to the biophysical characteristics of the forest 
canopy (Yong, Sun, & Li, 2004). Such measurements include total cano- 
py height, height to sub canopy layers, heights associated with different 
percentages of pulse energy return, height of median energy (HOME), 
canopy density (if assumptions are made concerning ground/canopy 
reflectivity ratio), and canopy height variability (Sun et al., 2008; Yong 
et al., 2004). These structurally related measurements, in turn, can be 
related to forest biophysical characteristics of interest such as basal 
area, timber volume, aboveground biomass, and C stocks (Lefsky, 
Harding, Parker, Acker, & Gower, 2002; Sun et al., 2008). 

2.4. ASTER GDEM data 

GLAS pulse interactions with vegetation on topography with a 
significant slope can result in pulse broadening, which confounds 
the interpretation of the influence of the vegetation alone on the 
waveform (Harding & Carabajal, 2005). We processed Advanced 
Spaceborne Thermal Emission and Reflection (ASTER) Global Digital 
Elevation Map (GDEM) Version 1 (VI) data to reduce the impact of 
pulse broadening. All poor quality GLAS pulses in North America, 
western Eurasia and eastern Eurasia were eliminated using techniques 
similar to those described by Lefsky, Keller, Pang, de Camargo, and 
Hunter (2007) using the thresholds provided in Table 2. 

The Ministry of Economy, Trade, and Industry (METI) of Japan and 
NASA jointly constructed the GDEM VI data. Multiple along-track stereo 
acquisitions were processed from the near-infrared spectral band and its 
nadir-viewing and backward-viewing telescopes with a base-to-height 
ratio of 0.6 at 15 m resolution. The GDEM data, dating back to 1999, 
covers 98% of the Earth's land surface in -1° tiles (60 x 60 km ground 
area). ASTER GDEM VI data were selected over Shuttle Radar Topogra- 
phy Mission (SRTM) data because of the higher resolution of the product 
(30 m vs. 90 m), and because of the near global coverage. SRTM data are 
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Table 2 

GLAS and ASTER GDEM metrics and valid ranges used to identify and exclude poor 
quality laser pulses. 


Metric 

Valid range 

AGB 

<500 

wflen 

0 < 80 

meanH, qmch, medH, Centroid, hi 4 , h 10 o 

-10< 50 

satndx 

<2 

ASTER GDEM slope 

<20° 


Variable definitions: 


AGB above ground live dry biomass, in megagrams per hectare 

wflen total vertical extent of the GLAS waveform from signal start to signal end 
meanH mean canopy height (Lefsky, Harding, Cohen, Parker, & Shugart, 1999) 
qmch quadratic mean canopy height profile (Lefsky et al., 1999) 
medH median of the canopy height (Lefsky et al., 1999) 
centroid height of the waveform centroids (Lefsky et al., 1999) 
hi 4 distance from signal beginning and ground peak, directly from GLA14 

h 10 o distance from signal beginning and ground peak, calculated from GLA01 

waveform 

satndx waveform saturation index, i.e. the number of within-canopy bins that 
have brightness values indicating the detector was saturated 
ASTER GDEM slope slope based on elevation A in a 3 x 3 GDEM window centered on 
the GLAS pulse. 


limited to areas below 60° N; the ASTER DEMs cover nearly all terrestrial 
surfaces (Hayakawa, Oguchi, & Lin, 2008; Hengl & Reuter, 201 1 ). GDEM 
VI accuracy is claimed to be 20 m at 95% confidence in the vertical di- 
mension and 30 m at 95% confidence in the horizontal dimension. For 
further reference, Toutin (2008) has provided detailed information 
about ASTER GDEM processing. 

Our processing developed mosaics of regional blocks of data 
(>500 tiles) using IDL/ENVI software. Blocks were processed with 
ENVI's topographic modeler, calculating min/max slope within a 
90 x 90 m window. Elevations <3 m and >6195 m were excluded 
from slope calculations. Limitations to using this data are analogous 
to those of SRTM as both report canopy elevations. 


2.5. Landsat and MODIS land products for stratification 

We used three different land cover maps to identify land cover 
types for stratification purposes within eco-regions of the boreal 
biome in North America and Eurasia. Land cover data for Canada 
was provided by the Earth Observations for Sustainable Development 
(EOSD) project, which generated a 25-m land cover map based on 
Landsat-7 Enhanced Thematic Mapper plus (ETM + ) data (centered 
on year 2000, with images from 1999 to 2002) released in 2006. 
Trees in this classification were defined as vegetation having a capac- 
ity to grow to heights > 5 m, and user accuracy for conifer forests was 
reported to be 86% in western Canada (Wulder et al., 2008). We have 
successfully used these data in a previous study of forest biomass in 
Quebec (Boudreau et al., 2008). 

In Alaska, we used the 30-m National Land Cover Data (NLCD) 
2001 classification based on Landsat-5 Thematic Mapper (TM) and 
Landsat-7 ETM + data collected between 1999 and 2004, with a ma- 
jority of the imagery between 1999 and 2002 (Selkowitz & Stehman, 
201 1 ). This provided a consistent, up-to-date land cover layer for the 
entire state. The NLCD land cover for Alaska is available through the 
Multi-Resolution Land Characteristics Consortium (MRLC). The 
MRLC purchased three dates of Landsat imagery for the entire U.S. 
and coordinated the production of a comprehensive NLCD land 
cover database (Homer et al., 2007). The MRLC consortium is specifi- 
cally designed to meet the current needs of federal agencies for na- 
tionally consistent satellite remote sensing and land cover data. 
Similar to EOSD, forests in this classification were defined as having 


heights > 5 m. Overall, accuracy of the product for forest non-forest 
is 83.9% (standard error of 2.1%) (Selkowitz & Stehman, 2011). 

For Eurasia, we adopted a similar approach for stratifying land 
cover. However, we used the International Geosphere Biosphere Pro- 
gram (IGBP) classification scheme from the 2004 500-m Moderate 
Resolution Imaging Spectroradiometer (MODIS) global land cover 
(MOD12Q1) data (Justice etal., 1998), instead of Landsat ETM + land 
cover maps to stratify cover types within sub-biome. The MOD12Q1 
data contains 17 land cover types; five of which are forest defined 
as having heights > 2 m (Friedl et al., 2010). Discrepancies in the 
identification of forested area near the taiga-tundra transition zone 
among remote sensing data products have been noted; however, 
good agreement exists throughout the boreal forest regions of north- 
ern Eurasia when comparing GLC-2000, GLOBCOVER, and MODIS 
IGBP classifications (Herold, Mayaux, Woodcock, Baccini, & 
Schmullius, 2008; Pflugmacher et al., 2011). Overall, global accuracy 
of the product has been stated to be 74.2% (Friedl et al„ 2010), and 
we used five of the forest classes plus woody savanna, which is a 
dominant cover type over much of far eastern Siberia and the higher 
latitudes of the North American boreal forest. 

We updated the two land cover products with recent burned area 
extents. We used a combination of fire polygon data obtained from 
the Canadian Forest Service and Alaska DNR and supplemented 
these with MODIS MCD45 lxl km burned area fire products as 
needed (Roy, Boschetti, & Justice, 2006). NRCan-CCRS-CFS in 
Canada have a program to use the MODIS burn mapping to prioritize 
areas for collection of Landsat (or SPOT) data to map the fires with 
greater spatial detail than that available from MODIS. MODIS data 
were used exclusively for the Eurasia portion of the boreal forest. 

2.6. Methods to develop biomass models 

2.6.1. Ecoregion stratification 

For a consistent biome-level analysis, we evaluated land area associ- 
ated with ecoregions contained within the boreal forest biome of North 
America and Eurasia. The boreal land area, and its component 
ecoregions, was defined by the World Wildlife Fund's ( WWF) ecoregion 
map of the world (Olson et al., 2001 ). We used these ecoregions, along 
with satellite-based land cover data, to create land cover strata, where- 
by the same land cover falling in different ecoregions was distinguished 
as unique strata. This stratum designation allowed for further refine- 
ment of land cover data and provided a means of grouping ground 
plots and LiDAR data (Fig. 1). Airborne measurements in Alaska and 
field measurements in northern Siberia enabled us to extend our ap- 
proach into these areas. 

2.6.2. Geographic data layers 

We have developed a data acquisition strategy to obtain our esti- 
mates of aboveground biomass and C for the boreal forest of North Amer- 
ica and Eurasia. To evaluate the variance of biomass estimates our basic 
study design relies on the following sources of data: 1) existing ground 
sample plots from various regions coupled with measurements of above- 
ground biomass on ground plots co-located with specific GLAS pulses; 
2) PALS and ALTM airborne LiDAR data acquired over ground sample 
plots, with additional airborne PALS and ALTM data acquired along se- 
lected GLAS orbits; 3) automated quality-screened L3c and L3f GLAS 
data for the boreal forests of North America — 930,400 pulses, automat- 
ed quality-screened L2a and L3a GLAS data in Scandinavia and western 
Eurasia — 738,159 pulses, and automated quality-screened L3c and L3f 
GLAS data in eastern Eurasia — 1 ,768,539 pulses (using criteria provided 
in Table 2); 4) Landsat land cover product for Canada (EOSD product), 
Alaska (NLCD product for the US; 30-m resolution) and MODIS 
MOD12Q1 land cover products for Eurasia; 5) ASTER GDEM VI (98% 
global coverage); and 6) WWF ecoregions. The resolutions and sources 
of these data are summarized in Table 3. 


278 


C.S.R. Neigh et al. / Remote Sensing of Environment 137 (2013) 274-287 


180' 





Table 3 

Data used to stratify study including source, resolution, and date. 


Dataset 

Source 

Res. (m) 

Date 

MODIS 

MOD12Q1 

(Eurasia) 

http://reverb.echo.nasa.gov/ 

500 

2004 

MODIS 

MCD45 

http://reverb.echo.nasa.gov/ 

1000 

2003-06 

ASTER 
GDEM VI 

http://reverb.echo.nasa.gov/ 

30 

2000-09 

Landsat MLRC 
(Alaska) 

http://www.mrlc.gov/ 

30 

2001 

Landsat EOSD 
(Canada) 

http://tree.pfc.forestry.ca/ 

25 

2000 

WWF 

ecoregions 

http://worldwildlife. 

org/publications/terrestrial- 

ecoregions-of-the-world 

o 

o 

o 

2000 


Minimum mapping unit applied after conversion from a vector ArcGIS shapefile. 


2) flying PALS along GLAS orbital tracks to measure profiling LiDAR 
height and canopy closure on individual GLAS pulses; and 3) developing 
a second set of models that relate PALS estimates of biomass to GLAS 
metrics. Once GLAS models for predicting forest biomass were calculat- 
ed for the different land cover strata, the satellite LiDAR can be used as a 
sampling tool to attribute land cover strata across the continent. 

2.8. Model based biomass estimation 

The models reported below are similar to those reported by Nelson et 
al. (2012) and Stahl et al. (2011). More information concerning the for- 
mulation of these models can be found in those references. In our design, 
we use the models below in a three-phase design in North America and 
western Eurasia (ground-airborne-GLAS with PALS and ALTM airborne 
LiDAR) and a two-phase design (ground-GLAS) in eastern Eurasia. An 
estimator of the mean value fi Yc (e.g., biomass, in Mg ha -1 ) in stratum 
c is: 


f*Y c = 


^E£iQc(«c) 


where G jc (d c ) is the total of all the per-hectare biomass estimates on 
GLAS shots that intercept stratum c along GLAS orbit i, where m is the 
number of GLAS orbits crossing that area of interest, and n ic equals the 
number of GLAS shots intercepting stratum c on orbit i. 

The mean per-hectare biomass across all C strata, can then be mul- 
tiplied by the known total area to obtain an overall total which is: 


fly — Ec=i W c fi Yc 


Boreal Forest 
Land Cover 


Forest 

□ Evergreen Needleleaf 
I Deciduous Needleleaf 

Deciduous Broadleaf 
Mixed 

Shrubland 

□ Closed 
I I Open 

Othe r Cover Types 

Woody Savanna 

□ Grassland / Savanna 
Permanent Wetland 

I Cropland 


PALS flight lines 
Boreal Forest 


Field & Airborne 


Sampling 


Number of 
field plotsl 


1,000 

150 

10 


Application of Carbon 


Equations 


Ecoregion groups 
North ^merica 

5 See Table 4 for 

Eurasia corres P on( ^ in Q 
^ equations 

2-Phase AGB model 
(Ground-GLAS) 


Fig. 1. Circumboreal forest stratification, a. Land cover information derived from 
MODIS and produced in the International Geosphere Biosphere Program classification 
scheme, b. The map shows the biomass sampling density throughout the study area 
with number of field plots shown by the size of gray circles. Red depicts portable air- 
borne laser system (PALS) flight lines over co-located Geoscience Laser Altimeter Sys- 
tem (GLAS) shots, c. A color-coded map depicting the regional spatial stratification of 
biomass models applied. 


where W c is the GIS-based area proportion of stratum c and 

£c=,W c =1.0. 

The variance of this total can be estimated as follows: 


, 7 _ 1 y- y- W c W d J2l^G ic (a c )-fi Yc n ic )(G id ( a d )-p Yd n id ) 
W ml sh n c n d m - 1 

« « W c W d |s.P!' ( — 

+ ££-=^££cov S2 ( 

c=ld=l n c n d j, h ' 


*hd, 


G jlc G j 2 d 


( 3 ) 


2.7. North America ground-PALS-GLAS sampling 


The three-phase design used in North America involved: 1 ) flying 
PALS over existing USFS-FIA ground plots and the Canadian plots to 
develop models to relate ground-measured biomass to PALS metrics; 


where n c is the denominator of model 1, i.e., the coverage number of 
GLAS shots that intercept stratum c on a given flight line. 
Cdv S2 (aj iC ,6tj 2d ) is the covariance of the ji and j 2 coefficients in the 
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predictive model used to predict biomass in strata c and d, and p c and 
p d are the number of parameters, i.e., the number of independent var- 
iables, intercepts in the models used to predict biomass for strata c 

m T 

and d. G' c = m X/ (*it> a ) where gj (x if , a) is the value of the 1st 

i t=i 

derivative with respect to the jth coefficient of the predictive biomass 
model for the tth GLAS shot in the cth stratum, ith orbit. 

The first term in model (3) describes sampling error; the second term 
(second line of model (3)) describes model error. Model error character- 
izes the variability of the coefficients of the predictive model(s) and de- 
scribes how those coefficients and hence the predictive line or surface 
would change with repeated sampling. The model error in (3) does not 
describe prediction error, or the variability of a predicted biomass value 
around one specific prediction model. Prediction error could be added, 
as per Nelson et al. (2012), but it was not considered in this investigation. 

Model (3) applies to situations where the models used to predict 
biomass in stratum c and d are: 1 ) dependent (e.g., different models 
developed using the same X, Y observations) and/or; 2) are one and 
the same, i.e., where the same model is used to estimate biomass in 
more than one cover class. If, as in the current study, models are de- 
veloped independently, then the second term in Eq. (3) simplifies to 
the following: 


- £ £ £ Cdv(a JiC , &j c ) Gj c G hc . 

c=i n c Jt j2 \ ' 


For an individual stratum, the variance estimator is: 


^ ^ £j=1 (Qc( a c) f*Yc n ic) | ' \- \ - Cnyf Qf a ) C' C 

m(m-l) + n 2 c ui^ C \ hc ' hc ) Cj ’ c GjjC ' 


( 4 ) 


3. Results 

3.1. AGB models 

To use GLAS as a sampling tool to estimate C, predictive models must 
be developed relating GLAS heights to ground measurements of AGB 
density. For the purposes of reporting results of total C mass and densi- 
ty, we divided our AGB model estimates by 50% (Houghton et al., 2000). 
In North America, an airborne LiDAR profiler was flown over existing 
ground plots where AGB had been measured as well as along 
-15,000 km of GLAS orbits from Alaska to Quebec. We developed two 
sets of models, the first relating ground-measured AGB to airborne 
LiDAR height and canopy density metrics, and the second set relating 
airborne LiDAR estimates of AGB to GLAS metrics — a three-phase de- 
sign. Western Eurasia followed a similar three-phase design, which 
employed an airborne LiDAR scanner in place of the profiler. GLAS foot- 
prints and field plots were overflown with ALTM on an area in southern 
Norway. The ground, airborne LiDAR, and GLAS measurements were 
used to formulate the models needed to generate biomass predictions 
for western Eurasia. Eastern Eurasia employed a two-phase approach 
relating field measurements directly to the GLAS measurements with- 
out the airborne LiDAR intermediary. 

AGB models were developed for five circumpolar regions and, 
within each region for selected ecozones, for up to five independent 
strata for North America including, wetlands, hardwood, conifer, 
mixedwood, and burned areas. Two independent strata in western 
Eurasia were identified — conifer and hardwoods. In eastern Eurasia, 
we applied one generic model to all five strata due to the limited 
number of field plots in each of the strata. To develop these predictive 
models, we attributed every co-located GLAS footprint with data 
layers. We analyzed all possible subsets of GLAS metrics to identify 


those independent GLAS variables that best predicted AGB estimates 
based on LiDAR profiler measurements in North America, LiDAR scan- 
ner measurements in western Eurasia, and field measurements in 
eastern Eurasia. Once AGB models were established, we applied 
them to GLAS pulses by region, ecozones, and strata to provide C es- 
timates across the entire boreal forest. Table 4 summarizes the 
models with variable definitions by region, ecozones, and strata that 
were used to generate the AGB estimates. 

3.2. AGB estimates 

We mapped by stratum the spatial distribution and standard error of 
C across the boreal biome, excluding shrubs (Fig. 2). Shrubs were ex- 
cluded for two reasons: 1) the lack of boreal shrub allometric models 
and plots located in shrub areas; and 2) the lack of sensitivity of GLAS 
to vegetation < 5 m in height (Nelson, 2010). We estimated total 
biome C mass to be 38.0 ±3.1 Pg (petagrams) distributed across a for- 
ested region encompassing 11.9 M km 2 . We report our C estimates by 
stratum in Table 5, within each of five general circumpolar regions 
including Alaska, western Canada, eastern Canada, western Eurasia 
(west of the Ural Mountains), and eastern Eurasia (east of and including 
the Urals). Mean C density for the entire boreal forest biome is 20.3 ± 
1.7Mgha~' (million grams per hectare), with eastern Eurasia 
reporting the largest per-hectare quantity, 22.7 ± 3.4 Mg ha -1 . The 
eastern Eurasia standard errors are large relative to the standard errors 
reported for the other four regions due to the approach we used to pre- 
dict biomass. A two-phase sampling approach was employed in eastern 
Eurasia, an approach that entails the use of a “direct” ground-GLAS 
model that predicts ground-based estimates of biomass as a function 
of GLAS measurements. In the other four regions, an airborne LiDAR 
was used as an intermediaiy to tie ground-based biomass estimates to 
GLAS metrics. The variance calculations in these three-phase designs 
ignore the ground-airborne LiDAR component and only incorporate 
the airborne LiDAR-GLAS model variances. The ground-airborne 
LiDAR model component was not included because our model-based 
estimator is limited to two-phase scenarios; we do not know how to 
legitimately incorporate the second (in this case, ground-airborne 
LiDAR) model uncertainty. We include in our variance estimate that 
model uncertainty which we believe is the larger of the two, i.e., the air- 
borne LiDAR-GLAS model. The within- and across-strata variances in 
the four regions that utilize a profiling or scanning LiDAR are therefore 
under-reported by an unknown amount. 

Across the entire circumboreal biome, 53% of the total C mass is in 
conifer forest (20.1 ± 2.4 Pg), 37% is in mixedwood forest (14.1 ± 
0.6 Pg), 5% is in hardwood forest (1.8 ± <0.1 Pg), 4% is in woody 
wetlands (1.6 ± 0.3 Pg), and 1% is in forest that burned between 
2000 and 2005 (0.4 ± 0.1 Pg). We note the burned area reported 
here is underestimated relative to other studies in North America 
(Kasischke et al., 2011), resulting in comparatively low C estimates 
for the burned area strata. 

3.3. How do these results compare with other C estimates? 

We performed validation of C estimates through study inter- 
comparison on a country basis for those that completely fall within 
the boreal biome. For Canada, we used data reported by the two 
Canadian inventory systems, CanFI and NF1, as a source of validation 
of our LiDAR-based C estimates. Within the managed forest area, we 
compared our estimates to those produced by the Carbon Budget 
Model of the Canadian Forest Sector (Kurz & Apps, 1999). In Scandinavia 
and Russia, validation was more problematic. For instance, Norway, prior 
to 2005, inventoried only its commercial forest and ignored significant 
stretches of mountain forest and northern forests that were not cur- 
rently marketable. Russian Federation forest inventory information 
tends to be in the southern boreal zone (similar to Canada), dated, 
and driven by local forestry needs. Our study and the Food and 


Table 4 

Study region model equations predicting aboveground live dry biomass (AGB). We approximate carbon to be 50% of AGB (Houghton et al., 2000). Models were developed via a three-phase sampling approach linking ground to airborne to 
satellite (Ground-PALS-GLAS) or a by a two-phase sampling approach linking ground to satellite (Ground-GLAS). Approaches varied by region (Fig. 1 ) in order to produce the most robust multiple linear predictive models with existing data. 
WWF ecoregions were grouped based on number of ground plots available, R 2 , RMSE, and VIFS <10, except for western Eurasia (D) stratum four where condition index was <12. Number of ground plots co-located to airborne measure- 
ments, and airborne to co-located GLAS acquisitions shown as (n). When constructed, condition number was the primary tool to quantify multicollinearity; post-analysis suggests that multicollinearity may be an issue here. For most of the 
boreal biome, models were developed for GLAS acquisitions L3c and L3f because regression models improved using these data as compared to using L2a and L3a. The exception was in western Eurasia, where L2a and L3a data were used. The 
Norway site was flown only over L2a. 


(A) Alaska WWF ecoregions 


Cook Inlet-Copper Plateau 3 

Interior Alaska-Yukon 

St. Elias and Coast 

Interior Mackenzie Delta 

Stratum 

Wetlands 

Hardwoods 

Conifer 

Mixedwood 

Burned b 


Ground-PALS model WWF ecoregion ID 

y p = 20.06 (h a ) - 1138 (h 70 ) - 4.12 50601, 50603, or 50604 

y p = 3.61 (hgs) - 6.94 (hgo) + 1.22 (deo) + 3.46 50607 

yl p = 34.15 (hgo) + 7.88 (h 10 ) - 5.51 51101 

y p = 12.28 (heo) + 17.05 (/ 130 ) + 17.08 511 1 1 or 51 116 

PALS-GLAS model 

% = 4.48 (h 14 ) - 627 (hso) - 2.10 (&*,) + 631 (acq3) + 11.66 
y g = 528 (hgo) - 434 (h 25 ) + 15.95 (acq3) + 37.16 
y g = 739 (hgo) - 8.75 (h 50 ) - 135 (acq3) + 36.98 
y g = 4.54 (h qc ) + 2.45 (h M ) + 18.96 (ocq3) + 5.76 
9g = -2.79 (h wflen _ ad] ) + 7.79 (hgo) - 8.05 (h 25 ) - 3.05 
(/slope) + 1622 (acq3) + 16.30 


RMSE 

Mg ha -1 R 2 

38.73 0.54 

39.15 0.58 

19.50 0.80 

29.24 0.74 

15.06 0.56 

23.75 0.53 

20.00 0.55 

29.62 0.54 

17.37 0.57 


VIFf n 

2.87 157 

1.85 111 

1.33 52 

4.12 41 

1.97 283 

1.85 176 

1.94 345 

4.92 156 

6.34 179 


(B) Western Canada WWF ecoregions 

Ground-PALS model 

WWF ecoregion ID 

RMSE 
Mg ha - 1 

R 2 

VIFf 

n 

Muskwa/Slave Lake 3 


50502, 50607, 50608, 50609, 50610, 
50612, 50613, 50614, 50617, 50802, 
51101, 51111, or 51116 





Conifer 

y P = 2.98 (hgo) + 7.02 (h 10 ) + 2.03 


32.40 

0.64 

1.19 

184 

Hardwood® 

9 P = 8.52 (h 50 ) - 0.50 (d 80 ) + 37.12 


30.18 

0.65 

1.74 

51 

Mixedwood 

ylp = 12.70 (h c ) - 139 (d 80 ) + 28.90 


35.23 

0.64 

1.70 

35 

All other cover types® 

y p = -6.23 (h a ) + 10.76 (h 10 ) - 2.06 


22.22 

0.84 

5.74 

36 


(C) Central and E. Canada WWF ecoregions 


Ground-PALS model 


WWF ecoregion ID 


Central Shield Forest 
Conifer 
Hardwood 
Mixedwood 
All other cover types® 
Eastern Transitional Forest 
Conifer 
Hardwood® 

Mixedwood 
All other cover types® 
Quebec and Maritimes 3 
Conifer 
Hardwood® 

Mixedwood 
All other cover type® 
Stratum 
Wetlands 
Hardwoods 
Conifer 
Mixedwood 
Burned d 


9 P = 12.04 (h 20 ) - 7.96 (s a ) - 1.18 
y P = 13.19 (hqc) - 0.73 (dm) + 2.36 
y P = 7.03 (hgs) + 0.22 (d4o) - 10.15 
y p = -6.23 (h a ) + 10.76 (h 10 ) - 2.06 

9 P = 10.68 (h g5 ) - 14.50 (s c ) + 3.83 
y P = 8.52 (hso) ~ 0.50 (d 80 ) + 37.12 
y p = 10.93 (ha) - 0.24 (d 30 ) + 0.32 
9 P = -6.23 (h a ) + 10.76 (h, 0 ) - 2.06 

ylp = 4.90 (hgo) + 6.30 (h 10 ) - 0.19 

y P = 8.52 (h 50 ) - 0.50 (d 80 ) + 37.12 

y>„ = 7.75 (heo) - 0.42 (d*>) + 43.65 

y p = -6.23 (h 0 ) + 10.76 (h 10 ) - 2.06 

Western Canada PALS-GLAS model 

9 g = 5.52 (h 14 ) - 5.17 (h 50 ) 3 2.22 (acq3) + 4.82 

y g = 233 (h wflen _ adj ) + 3.44 (h qc ) - 5.26 (acq3) + 35.68 

9 g = 5.31 (h 14 ) - 4.19 (h 25 ) + 4.03 

9 g = 4.11 (hqc) + 2.26 (h 14 ) - 2.50 (h 25 ) + 20.61 

y g = 5.46 (hqc) + 4.07 (ht Jr) + 22.77 (acq3) + 5.05 


50602 or 50616 


50406 


50605, 50606, or 50611 


50610, 50613, 50614, or 51116 


RMSE 

Mg ha -1 R 2 

22.21 0.80 

25.96 0.56 

26.74 0.60 

22.22 0.84 

37.05 0.56 

30.18 0.65 

28.63 0.48 

22.22 0.84 

27.90 0.69 

30.18 0.65 

30.82 0.49 

22.22 0.84 

22.23 0.52 

27.58 0.63 

21.28 0.68 

23.40 0.59 

27.15 0.66 


VIFf n 

2.15 60 

2.10 41 

1.02 98 

5.74 36 

2.24 29 

1.74 51 

1.46 82 

5.74 36 

1.75 66 

1.74 51 

2.19 88 

5.74 36 

1.94 1567 

3.12 779 

1.38 3042 

3.48 570 

1.47 31 
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Stratum 

Wetlands 

Hardwoods 

Conifer 

Mixedwood 

Burned d 

D) Western Eurasia WWF ecoregions 
Sarmatic Mixed Forest, Scandinavian and 
Russian Taiga, and Scandinavian Birch 3 
Stratum 

Conifer (mature spruce, med-high prod, sites) 
Conifer (mature pine, med-high prod, sites) 
Conifer (mature spruce/pine, poor prod, sites) 
Conifer (young spruce/pine, all site types) h 
Stratum 

Conifer, W. Wetlands, and W. Savanna 
Mixedwood and hardwoods 

(E) Eastern Eurasia WWF Ecoregions 
Central and Eastern Eurasia, and East Siberiaa 
Stratum 


All cover types 


Central/Eastern Canada PALS-GLAS model 

y g = 6.08 (h 14 ) - 1.29 (f slope ) - 1.20 

y g = 3.92 (h qc ) + 6.09 (h 90 ) - 5.44 (h 25 ) + 9.32 

y g = 2.27 (h 14 ) + 5.54 (/i 90 ) - 9.13 (h 25 ) + 2.40 

y g = 5.74 (h 14 ) - 0.92 (senergy) - 5.01 (h 25 ) + 14.63 

y g = 5.46 (h qc ) + 4.07 (/it Jr) + 22.76 (acq3) + 5.05 


Ground-ALTM model 8 

9. = 11.2024 (h/.oo) 1 ™' (dhS“'‘ 5 

9. = 7.4764 (cvl)°"“ (him) 0 - 6037 (dl„,)° s6s3 

9. = 1.6446 (h/so)" 359 {df, 0 ) , wu (W so)““ 529 (dfe,)-" 66 ' 

}, = 10.0093 (tl/go) 1 " 78 (d/ 0 ) 3JSSSS (i lfto)- 03014 (fi/ao)-' 33 ® 3 ® 

ALTM-GLAS model 

y w = 5.95 (h 90 ) - 5.00 (h 25 ) + 4.72 

y w - 8.84 (/i 7 s) - 5.09 (hjs) - 7.03 


Ground-GLAS model 


y e = 13.60 (/i 75 ) - 14.30 ( h 25 ) - 3.49 


50406, 50602, 50605, 50606, 50608, 9.93 

50609, 50802, 50407, 5041 0, or 50416 26.40 

23.56 
21.48 
27.15 


0.76 

0.79 

0.62 

0.70 

0.66 


1.36 

3.91 

4.66 

2.03 

1.47 


1348 

96 

850 

999 

31 


WWF ecoregion ID 


80436, 80608, or 81110 


RMSE (Obs. vs. Predicted) 

Mg ha -1 R 2 VIl* 

32.44 0.86 1.07 

18.22 0.88 7.21 

13.72 0.86 7.03 

15.51 0.91 25.32 


n 


51 

47 

51 

52 


33.61 0.46 2.65 212 

29.29 0.77 4.47 73 


WWF ecoregion ID 


80505, 80601, 80603, 80604, 80605, 
80606, 80607, 80609, 80610, 80611, 
80444, 80809, 80519, 80817, or 81105 


RMSE 
Mg ha -1 


58.47 


R 2 VII* 


0.60 1.80 


n 


55 


Variable definitions: 

• Dependent variables 


y p PALS estimates of total aboveground live dry biomass, in Mg ha -1 , a Ground-PALS linear model formulated by regressing ground estimates of biomass against PALS height and canopy density measurements. 

y g GLAS estimates of total aboveground live dry biomass, in Mg ha -1 , a PALS-GLAS linear model formulated by regressing PALS estimates of biomass against L3a and L3c GLAS measurements. 

j ) a ALTM estimates of total aboveground live dry biomass, in Mg ha -1 , Ground-ALTM log model linearized by taking the natural log of both sides of the model, formulated by regressing ln(biomass) against ALTM height and density 

measures and then back-transforming. 

y w GLAS estimates of total aboveground live dry biomass, in Mg ha - \ a ALTM-GLAS linear model formulated by regressing ALTM estimates of biomass against L2a and L3a GLAS measurements. 

y e GLAS estimates of total aboveground dry biomass, in Mg ha -1 , a ground-GLAS linear model formulated by regressing ground estimates of biomass against L3c and L3f GLAS measurements. 

• Independent variables 
- PALS variables in ground-PALS models 

h a , h c , hqc average height, all pulses; average height, canopy hits only (>1.3 m); quadratic mean height, canopy hits, respectively 
bio. h 20 , / 130 . /iso. I> 60 , / 170 . hgo, hg 5 10th, 20th, ... 25th, 90th, and 95th height percentiles of the first return canopy height distribution 

d 3 o, d 4 o, d 60 , dso the 60th percentile cumulative canopy density in percent. To calculate dgo on a plot or GLAS shot, (1 ) sort all pulses between 1.3 m and /i 9 s, (2) divide this sorted array into 10 equal height bins, and (3) calculate the per- 
centage of pulses included in an above a given height bin relative to the total number of pulses in the plot or GLAS shot. So do = (# pulses above 1.3 m) / (n„), where n a is the total number of pulses in the plot or GLAS shot. If 
hg 5 = 21.3 m so that each vertical bin is 2.0 m, then d 10 = percentage of pulses above 3.3 m (=1.3 + 2.0), and d 90 = the percentage of pulses above 19.3 m (=1.3 -1- 9 (2.0)) 
s a , s c standard deviation of h a ,h c , respectively. 

- ALTM variables in ground ALTM models 


hfso, li/go, h/ioo 80th, 90th, and 100th height percentiles of the first echo canopy height distribution 
/dso, hi 70 50th and 70th density percentiles of the last echo canopy height distribution 

d/o, d/io, d/ 4 o relative cumulative canopy height densities above the 0th (2 m), 1st, and 4th vertical height layers of the first echo height distribution, respectively 
d/ 50 , dl 70 relative cumulative canopy height densities above the 5th (2 m) and 7th vertical height layers of the last echo height distribution, respectively 

cvl coefficient of variation of the last echo height returns. 


(continued on next page) 
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hqc quadratic mean height of the waveform, calculated as the square root [X] (normalized amplitude in a given canopy height bin) x (height of bin) 2 ] 

h 2 5 , h 50 , h 75 , h 90 20th, 50th, 75th, and 90th percentiles of the waveform heights; e.g., h 75 is that distance from the height bin below which 75% of the total waveform energy resides to the center of the ground peak 

h u top of tree height, signal start to ground peak, GLAM product 

hw/ien height of the waveform, signal start to signal end, GLAM product 

fsiope front slope to surface energy ratio from waveform as calculated by Sun et al. (2008) 
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Agriculture Organization-Forest Resources Assessment (FAO-FRA) report 
total aboveground Pg of C for 2005 by country. Note that a majority of 
studies do not report aboveground C separately. Furthermore, we are un- 
aware of any other biome-based C maps with which comparisons can be 
made. To facilitate comparison with other studies whose maps and esti- 
mates of C stocks are reported for broad administrative units, we have ag- 
gregated boreal C results to the country level. 

The countries by which we summarize the boreal C totals are 
those countries that contain ecoregions that are part of the boreal 
biome. It is important to note, when comparing total C that there is 
a large variation in each study's estimate of forested area. We provide 
Table 6, comparing our estimates by country with those from other 
studies. We found the following differences between our estimates 
and FAO-FRA: Russian Federation 6.6%; Canada 28%; Estonia 50%; 
Latvia 50%; Finland <0.1%; Norway <0.1% and Sweden 10%. In all com- 
parisons, GLAS-based estimates were less than or equal to FAO-FRA 
estimates. 



4. Discussion 

Our approach utilizes model-based estimators and a two- or 
three-phase statistical sampling procedure (Nelson et al., 2012; 
Stahl et al., 201 1 ) to estimate C in forests. This approach is the culmi- 
nation of long-term research that has focused on accounting for, and 
when possible, mitigating a variety of sources of uncertainty. Sources 
of uncertainty, or variance, include sampling variability, model error 
(i.e., variability of the coefficients), and the covariability among stra- 
ta across all GLAS orbits. An early airborne LiDAR sampling study in 
Delaware (Nelson, Short, & Valenti, 2004) used Line Intercept Sam- 
pling (LIS) techniques (Kaiser, 1983) to calculate land cover stratum 
estimates and variances. LIS assumes that the LiDAR flight lines are 
randomly allocated across the study site when in fact, they were sys- 
tematically spaced. Subsequent work in Quebec (Nelson, Boudreau, 
et al., 2009) and Siberia (Nelson, Ranson, et al., 2009) (1) accounted 
for regression error, (2) took into account the co-variability between 
adjacent and near-adjacent flight lines, and (3) tried to mitigate the 
potential inflationary effect of treating a systematic sample as a ran- 
dom sample by introducing successive-difference variance estima- 
tors. In Hedmark County, Norway we (Ntesset et al„ 2011) and 
others (Gregoire et al., 2011; Stahl et al., 2011) have further refined 
variance estimators in the context of sampling large areas with air- 
borne LiDARs. Gregoire et al. (2011) report on a model-assisted 
design-based statistical approach which may be employed in situa- 
tions where a sample of airborne LiDAR flight lines is collected in 
association with a probability-based ground sample, e.g., an array 
of national forest inventory plots. The technique set forth by Stahl 
et al. (2011 ) does not require the presence or integration of a proba- 
bility ground sample, making this model-based technique more 
amenable to assessment of remote, inaccessible areas that lack a 
design-based sample of ground plots such as the circumboreal 
areas targeted in this study. Eastern Eurasia had only one ground- 
GLAS model primarily because of limited field sampling, and in part 
because of the absence of any available airborne LiDAR. The large 
uncertainty in this region is evident in Fig. 2 with eastern Eurasia 
containing error that is excluded in the three-phase design. Errors 
for this region are two to three times as much as reported in the 
three-phase design. Future studies should integrate airborne with 
additional field data in remote areas of Siberia to reduce this error 
and develop a solution to account for model error in situations 
where two models are used in a three-phase design to estimate 
biomass. 

Our circumboreal assessment combined multiple data products 
assembled at different points of time and across many scales. This ag- 
gregation of data introduces many sources of error, some of which we 
do not account for, and these errors may propagate in unknown ways. 
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Fig. 2. Circumboreal forest, a, Aboveground carbon, b, c Error estimates of carbon den- 
sity. The images illustrate stratum-level estimates reported in Table 5. 


Sources of error that we believe could be better controlled, mitigated, 
or accounted for in future studies include: 

1) temporal error associated to non-coincident ground, airborne, 
and GLAS measurements by land cover stratum — expected to 


be small in boreal forests, and larger for warmer, faster growing 
forests; 

2) allometric error in models used to develop field plot biomass 
estimates — expected to be small in regional/national surveys 
(McRoberts & Westfall, accepted for publication); 

3) geolocation error between ground plots and airborne LiDAR (5- 
15m) — effect expected to be large as geolocation error exceeds 
10 m (Gobakken & Naesset, 2009); 

4) spaceborne LiDAR geolocation error (-15 m) — expected to be 
large; 

5) land cover product geolocation error coupled with minimum 
mapping units, differences in forest height definition, and crown 
closure thresholds — unknown; 

6) digital terrain model artifacts impacting slope calculations and 
geolocation errors — unknown; and 

7) land cover product classification error. These sources of error we 
believe affected our results by an unknown amount. Future stud- 
ies should consider these sources of error to improve C estimates 
— expected only to minimally result in regional across strata inac- 
curacy but it may result in assignment of biomass to incorrect 
strata. 

These sources of error, we believe, affected our results by an un- 
known amount. Future studies should consider mitigating or account- 
ing for these sources of error to improve C estimates. 

Our study design, developed from prior investigations, imposed a 
sampling structure that we believe maximized utility of the products 
available to estimate circumboreal C. Using GLAS as our sampling 
tool, we characterize strata defined by publically available land 
cover maps. In Alaska, that was the 2001 30-m NLCD product that 
defines trees or forest as being vegetation that attain heights > 5 m 
(Gobakken & Ntesset, 2009; Homer et al., 2007). In Canada, we 
employed the circa 2000 25-m EOSD product where forests were 
defined as vegetation with a capacity to grow > 5 m (Wulder & 
Nelson, 2003). In Eurasia, lacking wall-to-wall ETM + land cover 
maps, we employed the 2004 500-m MODIS MCD12Q1 IGBP product 
where forests were defined as vegetation > 2 m (Friedl et al., 2002). 
Future studies in Russia may find the hybrid land cover dataset useful 
as it provides more cover classes at higher resolution (Schepaschenko 
et al., 2011). 

These land cover maps sources provide a stratification framework 
but do not limit the vegetation heights below which we cannot mea- 
sure tree heights or estimate biomass with GLAS. Our vegetation 
height sensitivity is not driven by map definition; it is actually driven 
by the definition of tree used by the ground plot surveys. Typically, 
ground crews use minimum dbh and height limits to define what 
constitutes a tree. These dbh limits are approximately 5-10 cm, and 
1.3 min height in Norway for example (Stahl et al., 201 1 ). Depending 
upon the particular ground sample being discussed these minimum 
thresholds vary. GLAS is relatively insensitive to short sparse forest, 
but within our forested map strata as defined by the NLCD, EOSD, 
and MCD12Q1, we use GLAS nonetheless to estimate biomass across 
the entire height and density gradient. We purposely excluded 
non-forest strata, e.g. shrubs, water, shrub wetlands, that may in 
fact support some forest biomass and assume that if Landsat or 
MODIS identifies a particular GLAS shot as falling in the shrub stra- 
tum, then the biomass assigned to that GLAS shot is equivalent to 
zero, regardless of whatever vegetation height might have been mea- 
sured by GLAS. We do this based on our experiences in Quebec, where 
GLAS “measured" significant biomass in areas that we knew were ex- 
tensive areas of relatively flat rock, i.e. the Precambrian Shield in 
northern Quebec (Boudreau et al., 2008; Nelson, Boudreau, et al., 
2009). By applying this approach, we believe that we probably miss 
significant amounts of forest biomass in areas defined by land cover 
maps as non-forest, but GLAS is not the appropriate tool to infer bio- 
mass in these areas. 
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Table 5 

Aboveground C of the boreal forest, including sub-regions, and strata. 



North America 



Eurasia 




Alaska 

W. Canada 

E. Canada 

West 

East 

Circumboreal 

C by region 

Average density, Mg ha -1 

14.0 

15.8 

20.2 

19.4 

22.7 

20.3 

Standard error 

0.3 

0.3 

0.8 

1.0 

3.4 

1.7 

Total biomass Pg 

1.1 

5.1 

4.8 

6.3 

20.8 

38.0 

Standard error 

<0.1 

0.1 

0.2 

0.3 

3.1 

3.1 

Area 

Forest area, km 2 x 10 3 

370.0 

1779.6 

1547.1 

2350.4 

5881.2 

Total area 
11,928.3 

Non-forest area km 2 x 10 3 

382.2 

1433.3 

820.4 

891.3 

3278.4 

6805.6 

Percent forest 

49.2 

55.4 

65.4 

72.5 

64.2 

63.7 

Percent of boreal forest 

4.0 

17.2 

12.6 

17.3 

48.9 

100.0 

Area by stratum, km 2 x 1 0 3 

Wetlands 

58.7 

480.7 

300.2 

79.0 

331.4 

1250.0 

Hardwood 

49.5 

182.9 

91.3 

17.8 

28.5 

370.0 

Conifer 

205.9 

921.2 

780.9 

1363.6 

4011.7 

7283.3 

Mixedwood 

45.8 

163.7 

365.0 

887.4 

1379.1 

2841.0 

Burned 

10.2 

31.1 

9.7 

2.5 

130.6 

184.1 

C by stratum (Tg) 
Wetlands 

69.3 

822.4 

354.2 

101.0 

264.7 

Total C 
1611.6 

Standard error 

6.2 

22.8 

19.4 

21.3 

256.6 

259.3 

Hardwood 

216.2 

822.8 

518.1 

53.0 

172.4 

1782.5 

Standard error 

6.5 

16.7 

18.6 

5.0 

13.7 

29.6 

Conifer 

579.6 

2707.3 

2310.9 

3499.0 

11,053.2 

20,150.0 

Standard error 

13.7 

64.6 

162.9 

251.4 

2416.1 

2435.5 

Mixedwood 

173.2 

657.0 

1568.0 

2623.3 

9041.1 

14,062.6 

Standard error 

10.1 

15.2 

60.8 

208.5 

582.9 

622.3 

Burned 

16.6 

76.2 

26.1 

4.1 

227.2 

350.2 

Standard error 

1.3 

9.8 

4.8 

1.4 

89.9 

90.6 


5. Conclusions 

Our circumpolar estimates of C, with standard errors, within the 
boreal forest build on the approaches above, and thus represent a sig- 
nificant advance in terrestrial C accounting. Furthermore, the sam- 
pling approaches and statistical analyses that we used provide 
important guidance to design and use of future spaceborne LiDAR 
systems for statistically evaluating terrestrial aboveground C stocks. 
The ICESat-2 photon-counting LiDAR, due for launch in 2016, and 
the recently approved BIOMASS satellite mission using P-band 
(435 MHz) SAR (Le Toan et al., 2011) could aid in re-measuring the 
circumboreal forest. (We note that the ability to collect P-band 
radar data over Europe and North America is currently not permitted 
at this frequency due to military restrictions; de Selding, 2013. The 
launch of BIOMASS is not planned until 2020, with the topic currently 
under high-level discussion — the outcome of which will affect boreal 
C mapping interests and initiatives). In general, these new systems 
and measures are desired to decrease the uncertainty around current 
estimates and to aid in assessing whether C management activities 
aiming to increase C stocks have indeed had a measurable 
large-area impact. Spaceborne LiDAR also has the potential to serve 
as a monitoring tool for quantitatively assessing how aboveground 
C stocks will vary with the large changes in climate predicted for 
northern regions (Soja et al., 2007). 

Future monitoring of boreal C stocks will evolve with next- 
generation, space-based LiDAR instruments to produce a time series of 
C inventories. These time-series will provide insight to the spatial and 
temporal distribution of the changes in forest C stocks, and will enhance 
our ability to model and monitor the C cycle. For present modeling 
studies, the numbers reported here, with error bars, will place an im- 
portant constraint on aboveground C pools that historically have been 
an under-constrained modeling problem. Understanding C changes 
resulting from forest management activities, climate change, or 


modification in the rates or intensity of harvest is critical to meet United 
Nations climate change treaty reporting goals. 
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Table 6 

Country comparison of boreal forest carbon (C) estimates. We report aboveground C ± 1 standard error. 


Country Forest area 

Total C 

Standard error 

Mean C density 

Standard error 

Year 

Source 

(10 6 ha -1 ) 

<pg) 

(pg) 

(Mg ha~') 

(Mg ha~') 



Russian Federation 







771.1 

23.6 

N/A 

N/A 

N/A 

1990 

Krankina and Dixon (1994) 

771.1 

35.1 

N/A 

41.1 

N/A 

1990 

Isaev, Korovin, Zamolodchikov, 
Utkin, and Pryaznikov (1995) 

886.6 

22.1 

N/A 

N/A 

N/A 

1993 

Krankina, Harmon, and Winjum (1996) 

882.0 

28.0 

N/A 

36.3 

N/A 

1998 

Alexeyev and Birdsey (1998) 

886.5 

26.5 C 

N/A 

38.6 

N/A 

1993 

Shvidenko and Nilsson (2003) 

523.6 a 

23.1 

N/A 

44.1 

N/A 

2000 

Houghton et al. (2007) 

826.6 a 

33.3 

N/A 

40.3 

N/A 

2000 

Houghton et al. (2007) 

N/A 

33.7 d 

N/A 

N/A 

N/A 

1990 

Goodale et al. (2002) 

642.2 

24.4 

N/A 

38.0 

N/A 

1999 

Dong et al. (2003) 

821.8 

36.0 

N/A 

N/A 

N/A 

2000 

Pan et al. (2011) 

845.6 

37.5 

N/A 

N/A 

N/A 

2007 

Pan et al. (2011) 

808.8 

25.8 C 

N/A 

N/A 

N/A 

2005 

FAO-FRA e 

725.2 

24.1 

±3.1 

33.2 

4.2 

2005 

This study f 

Canada 







229.7 

14.4 

N/A 

N/A 

N/A 

2000 

Pan et al. (2011) 

229.4 

14.0 

N/A 

N/A 

N/A 

2007 

Pan et al. (2011) 

404.2 

N/A 

N/A 

28.4 C 

N/A 

1989 

Kurz and Apps (1999) 

N/A 

12.9 d 

N/A 

N/A 

N/A 

1990 

Goodale et al. (2002) 

512.6 s 

9.7 s 

±2.0 S 

19.0 s 

4.0 

1990 

Botkin and Simpson (1990) 

239.5 

10.6 

N/A 

44.1 

N/A 

1999 

Dong et al. (2003) 

310.1 

13.7 C 

N/A 

N/A 

N/A 

2005 

FAO-FRA e 

332.7 

9.9 

±0.2 

29.7 b 

0.6 b 

2005 

This study 

Estonia 







2.3 

0.1 

N/A 

47.5 

N/A 

1999 

Dong et al. (2003) 

2.3 

0.2 C 

N/A 

N/A 

N/A 

2005 

FAO-FRA e 

2.3 

<0.1 

<±0.1 

23.3 

3.0 

2005 

This study 

Latvia 







3.5 

0.2 

N/A 

49.8 

N/A 

1999 

Dong et al. (2003) 

3.3 

0.2 C 

N/A 

N/A 

N/A 

2005 

FAO-FRA e 

2.8 

0.1 

<±0.1 

41.8 

3.6 

2005 

This study 

Finland 







17.2 

0.6 

N/A 

34.9 

N/A 

1999 

Dong et al. (2003) 

22.2 

0.7 C 

N/A 

N/A 

N/A 

2005 

FAO-FRA e 

27.9 

0.7 

<±0.1 

22.7 

2.0 

2005 

This study 

Norway 







7.0 

0.3 

N/A 

37.2 

N/A 

1999 

Dong et al. (2003) 

9.7 

0.3 C 

N/A 

N/A 

N/A 

2005 

FAO-FRA e 

10.7 

0.3 

<±0.1 

29.8 

1.6 

2005 

This study 

Sweden 







26.5 

1.1 

N/A 

39.9 

N/A 

1999 

Dong et al. (2003) 

28.2 

1.0 C 

N/A 

N/A 

N/A 

2005 

FAO-FRA e 

33.1 

0.9 

±0.1 

27.3 

1.7 

2005 

This study 


a Estimate of forest area from MODIS MOD12Q1 and GLC2000 land cover datasets. 

b C and standard error estimated as mean between East and West Canada World Wildlife Fund (WWF) ecozones. 
c C reported as living aboveground only biomass. 
d Live vegetation C reported with understory. 

e Food and Agriculture Organization Forest Resources Assessment (FAO-FRA) derived from 2005 database available online at: countrystat.org. 
f Our Russian Federation forest area estimate is boreal only. We did not include non-boreal forests in the southeast or southwest. 
g Values reported for all of boreal North America including Alaska. 
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